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ABSTRACT 

We show that measured velocity dispersions of dwarf spheroidal galaxies from about 4 to 10 km/s are 
unlikely to be inflated by more than 20% due to the orbital motion of binary stars, and demonstrate 
that the intrinsic velocity dispersions can be determined to within a few percent accuracy using 
two-epoch observations with 1-2 years as the optimal time interval. The crucial observable is the 
threshold fraction — the fraction of stars that show velocity changes larger than a given threshold 
between measurements. The threshold fraction is tightly correlated with the dispersion introduced 
by binaries, independent of the underlying binary fraction and distribution of orbital parameters. We 
outline a simple procedure to correct the velocity dispersion to within a few percent accuracy by using 
the threshold fraction and provide fitting functions for this method. We also develop a methodology for 
constraining properties of binary populations from both single- and two-epoch velocity measurements 
by including the binary velocity distribution in a Bayesian analysis. 
Subject headings: binary stars, theory — galaxies: kinematics and dynamics 



1. INTRODUCTION 

In recent years, a large number of dwarf spheroidal 
satellite galaxies of the Milky Way have been discovered 
using the Sloan Di gital Sky Survey (|Willman et al.ll2005l 
iZucker et~a l. 2006). These galaxies are much fainter than 
previously known Milky Way satellites, having larger 
mass-to-light ratios and velocity dispersions that rang e 
from 7.6 km/s down to 3.3 km/s (|Simon fc Gehall2007l ). 
Estimating the amount of dark matter contained in these 
galaxies becomes more susceptible to error than in the 
larger dwarf spheroidals. This is not only because of the 
statistical error associated with small stellar samples, but 
also because sources of velocity contamination constitute 
a larger fractional error due to the galaxy's small intrin- 
sic velocity dispersion. Potential sources of contamina- 
tion come from foregroun d Milky Way stars , atmospheric 
jitter in red giant stars (jPrvor et al.l 119881 ). and an in- 
flated veloci ty dispersion due to t he orbital motion of 
binary stars ((Olszewski et al.lll996| ). Among these, bina- 
ries have been the most difficult to correct for because the 
binary distribution of velocities in environments beyond 
the solar neighborhood is not well known and difficult to 
constrain without a large number of high-precision radial 
velocity measurements. 

The most prominent signature of binary stars in these 
galaxies is a high- velocity "tail" in the velocity distribu- 
tion due to short-period binaries, which gives the distri- 
bution a higher kurtosis than that expected from a Gaus- 
sian. Because of this departure from Gaussianity, the 
intrinsic velocity dispersion is usually estim ated by us- 
ing a robust estimator such as the biweight (| Beers et al.l 
1990), by discarding velocity outlier stars from the data 
sample, or by a combination of both approaches (cf. 
iMateo et al.lll99ll ) . While these techniques eliminate the 
largest component of the binary dispersion, binaries in- 



flate the first-order (Gaussian) c omponent of the disper - 
sion as well. In previous stu dies, lOlszewski et al.1 (|1996l ) 
and lHargreaves et al.1 (|1996l ) used Monte Carlo simula- 
tions to show that the dispersion estimated by these tech- 
niques is inflated due to binaries by an amount which is 
small compared to the statistical error in their data sets. 
Since that time, the number of stars in then-known dwarf 
spheroidals with spectroscopic data has i ncreased consid- 
erably , from less than 100 at the time of lOlszewski et afl 
(1996) to more than 10 00 for Draco, Fornax, and Ca- 
rina ([Walker et al. 2009]). For such large samples, we 
will show that the bias in the first-order dispersion due 
to binaries is larger than the statistical error. More im- 
portantly, for galaxies whose intrinsic dispersion is small, 
the first-order dispersion may be inflated by somewhat 
more than the statisical error even in data sets as small 
as 100 stars. 

To improve previous estimates of the intrinsic disper- 
sion for these cases, it is necessary to model the velocity 
distribution of the binary population and investigate the 
behavior of the binary dispersion as model parameters 
such as binary fraction are varied. With this approach 
it is desirable to find the best possible constraints on 
the binary velocity distribution, and observations taken 
at multiple epochs are very useful for this purpose. In 
addition, estimating the binary fraction and other bi- 
nary properties in galaxies and clusters beyond the so- 
lar neighborhood is useful in its own right, since this 
has been diffi c ult to predict accura tely in simulations 
(jTohlind [20021 iGoodwin et "ail l200l . Binaries also af- 
fect higher-order moments of the velocity distribution, 
e.g. kurtosis, which has been shown in principle to pro- 
vide useful informa tion about the mass profile of galaxies 
([Lokas et al.l 120051 ). 

In section [5] we derive an analytic formula for the 
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center-of-mass velocity distribution of binary stars. In 
section [3] we use this distribution to derive likelihoods 
for both single- and two-epoch velocity data. We demon- 
strate in section |4] how properties of a binary population 
can be constrained from multi-epoch data. Fig. |4] shows 
that repeat measurements with a baseline of 1-2 years 
is sufficient to constrain the binary population. We will 
prove that over a range of velocities of a few km/s, the 
binary fraction is nearly degenerate with the parame- 
ters characterizing the distribution of periods, and in sec- 
tion [5] we use this fact to develop a model-independent 
method for correcting the dispersion due to binaries by 
using multi-epoch data. Fig. |TD] encapsulates our main 
result, from which we conclude that the velocity disper- 
sions of dwarf spheroidals are unlikely to be inflated by 
more than 20% by binaries. The procedure to correct 
the dispersion for binaries is summarized at the end of 
section [5] In section |6] we outline a method to combine 
single- and multi-epoch data in a Bayesian analysis, and 
discuss the issue of foreground contamination and how 
this affects the apparent binary population. In section 
[7| we find a fitting function for the velocity distribution 
of binary stars and use it to derive an analytic formula 
for the binary dispersion. Finally, in section |8] we further 
explore the degeneracy of model parameters and discuss 
prospects for constraining the period distribution of bi- 
nary systems from multi-epoch velocity data. 

2. DISTRIBUTION OF VELOCITIES IN THE 
CENTER-OF-MASS FRAME OF BINARY SYSTEMS 

Spectroscopic velocity measurements of binary stars 
are typically dominated by the primary star due to 
the difference in luminosities and spectra between the 
two stars. This is especially true for stars that lie 
on the red giant branch, for which the luminosity- 
mass re lation steepens drast i cally. By way of com- 
parison, iDuauennov fe Mavorl (|199lD found the ratio of 
secondary-to-primary masses of binary stars within the 
solar neighborhood to peak at m/M ps 0.23 with only 
ps 20% of the stars having mass ratios larger than 0.6. 
For systems with a mass ratio of 0.6 and the primary ly- 
ing near the base of the red giant branch, the luminosity 
ratio would be at most I ps (m/M) 2 ' 3 ps 0.3; this ratio be- 
comes much smaller if the primary star is further along on 
the red giant branch. Thus, although the cross correla- 
tion function of the st ellar spectra may be double-peaked 
(jTonry fe Davisiri979[ ). the spectrum of the primary star 
will likely dominate the signal unless the mass ratio is 
close to 1, in which case the velocities of the two stars 
will be nearly equal. In view of the preceding arguments, 
for the remainder of this paper we will use the velocity of 
the primary star in modeling the spectroscopic velocities 
of binary systems. 

To model the velocity distribution of the primary star 
in a population of binary systems, first we must find its 
velocity distribution in the center-of-mass frame of the 
binary. The motion of two stars orbiting each other can 
be simply expressed in terms of four parameters: the 
semi-major axis a, eccentricity e, period P, and orbital 
angle <j> with respect to the center of mass. If we also 
specify the orientation of the orbital plane in terms of 
Euler angles, with the z-axis pointing along the line of 
sight, the line-of-sight velocity of the primary star in the 
center-of-mass frame is given by 
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g e (0, tp, <f)) = sin 9 [cos(V> — <f>) + ecos(0)] 
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where a is the semi-major axis of the primary star's or- 
bit, and 9 and ip are the second and third Euler angles, 
respectively. The azimuthal Euler angle does not appear 
explicitly because v z is invariant under rotations about 
the z-axis. 

Now taking the log of eq. Q] and subsituting Kepler's 
third law, we find 



logP = k - 3\og\v z 
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where P is in years, m is the mass of the primary star in 
solar masses, q is the ratio of secondary-to-primary mass, 
and v z is in units of AU/year. 

As an aside, we note that eqs. [3] and [4] can be used 
to make a back-of-the-envelope estimate of the velocity 
scale associated with a given orbital period. By averag- 
ing over orientations, mass ratio and eccentricity (whose 
PDF's are given below), we find 
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where M is the mass of the primary star. In dwarf 
spheroidal galaxies and globular clusters where red giants 
have masses M ps 0.8M Q , the above estimate shows that 
periods longer than a few decades will yield velocities less 
than 2 km/s. Meanwhile, velocities larger than 10 km/s 
will be dominated by binaries with periods shorter than 
1 month. The above estimate is somewhat sensitive to 
the mass ratio; whereas we used the approximate mean 
value q ps 0.4 in deriving eq. if a system has a mass 
ratio q ps 0.8, the coefficient in front becomes ss 10 km/s. 

To find a distribution in the center-of-mass velocity 
v z , the distribution of orbital periods must be averaged 
over all the parameters in eq. |4j We must also average 
over the time taken to traverse one orbital cycle, with 
all times being weighted equally. We express the orbital 
angle <$> in terms of the eccentric anomaly parameter 77, 
so that a uniform distribution in time corresponds to a 
distribution f(rj) = (1 — ecosrj)/2n. 

The orbital periods of G-dwarf stars in the solar nei gh- 
borhood were found by IDuauennov fc Mavorl (|1991D to 
follow a log -normal distributi o n wit h a mean period of 
180 years. iFischer fc Marcvl (|1992t) also found a log- 
normal period distribution for M-dwarfs in the solar 
neighborhood with a mean period similar to that of the 
G-dwarfs. In terms o f logarithm of the pe ri od P , the 
distribution found by IDuauennov fc Mavorl (| 1 9 9 If) has 
mean /ii og p — 2.23 and dispersion u\ og p = 2.3 (where P 
is in years, and the logarithm is base 10). We shall use 
this as the fiducial binary model in this paper, but will 
also allow /ii og p and <J\ og p to take on other values. 

As with the period distribution, we use distributions of 
the mass ratio q and eccentricity e observed in G-dwarf 
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stars in the solar neighborhood ([Duquennov fc Mavorl 
[1991 . When comparing different empiric ally derived 
mass functions, iDuquennov fc Mayor! (|199lD snowed the 
distribution of mass ratios q was best fit by the Gaus- 
sian mass function considered in iKroupa et al.l ()1990f ) 
with mean q — 0.23 and dispersion a q = 0.32. This 
is somewhat misleading, since for q > 0.5 the distribu- 
tion is in fact consistent with the power-law initial mass 
funct ions of Salpeter and Kroupa (Salpetcr 1955, Kroupal 
120011 ). At smaller mass ratios the distribution decreases 
sharply compared to the well-established Kroupa initial 
mass function of single stars, which implies that small 
mass ratios are strongly affected by interaction between 
primary and secondary during the formation of the bi- 
nary system. 

To further complicate matters, Ma zeh et al.1 ()1992[ ) 
found that short-period bin aries tend to have higher mas s 
ratios than those found bv IDuquennov fc Mavorl (|1991h . 
By analyzing spectroscopic binaries, they found that bi- 
naries with periods shorter than 3000 days have mass ra- 
tios consistent with a uniform distribution, although the 
Poisson errors i n their analysis were quite large. Sub- 
seque nt studies ([Goldberg et al. 2003, Halbwach s et al.l 
120041 ) have shown that the mass ratio distribution in 
short-period binaries is bim odal (also seen in previous 
samples; cf. iTrimbld I1990D . with a peak at low mass 
ratios similar to that of the long-period binaries but 
with an other peak near q 1. In the sample ana- 
lyzed by iGoldberg et al.l (|2003f ) , the peak at high mass 
ratios is smaller for primaries with masses larger than 
0.6 Mq, but larger for halo stars. In view of lingering 
uncertainties in the nature of the mass ratio distribu- 
tion, for si mplicity we will adopt the Gaussian distri- 
bution from IDuquennov fc Mavon (|1991| ) for long-period 
binaries (P > 3000 days) and a uniform distribution for 
short-period binaries (P < 3000 days). For the distribu- 
tion of the primary mass f(m) we will use the Kroupa in i- 
tial mass function corrected for binaries (|Kroupai r2002'). 

The distribution of ecc e ntricit ies /(e| logP) was found 
by IDuquennov fc Mavorl (|1991| ) to have three different 
regimes depending on the period. For periods of 11 days 
or shorter, the orbits are circularized due to tidal forces 
and are therefore approximated to have e — 0. For pe- 
riods between 11 days and 1000 days, /(e) can be ap- 
proximated by a Gaussian with mean e = 0.25 and dis- 
persion a e = 0.12. For periods longer than 1000 days, 
higher eccentricities are more prevalent and the distribu- 
tion approximately follows /(e) = fe 1 / 2 . Among these, 
the Gaussian regime (11 days < P < 1000 days) has 
the greatest impact on the velocity distribution at ve- 
locities of order km/s. While the adopted distributions 
of mass ratios and eccentricities are undoubtedly only a 
rough approximation to the true distributions, our cen- 
tral results (presented in section [5]) will prove to be quite 
insensitive to the nature of the adopted distributions. 

An important effect that must be taken into account 
is the effect of mass transfer between the stars if the pri- 
mary star is a red giant whose size is lar ger than the 
radius of its Roche lobe (jPacz vhskil 1 1 9 7 1| ) . In such a 
case, matter from the surface of the giant will accrete 
onto the smaller star and the separation between the 
stars will decrease. The end result is that either the 
other star will explode in a supernova la, or the stars will 



eventually m erge. This effect is no t inclu ded in the dis- 
tributions of IDuquennov fc Mavorl ()199lD because their 
sample consisted entirely of dwarf stars. Therefore, we 
make the approximation of excluding systems whose pri- 
mary star is larger than its Roche lobe, assuming the 
binary to be destroyed over a timescale much less than 
1 Gyr. We use an approxi mation to the ra dius r£,(a, q) 
of the Roche lobe given bv lEggletonl (fl98l . While this 
radius is not exactly correct for eccentric orbits, a re- 
cent smoothed-particle hydrodynamics si mulation of an 
eccent ric binary of mass ratio q = 0.6 by IChurch et al.l 
(2009) found the Roche lobe radius to decrease only 
slightly with eccentricity. They also derive a fitting func- 
tion for the Roche lobe radius with an eccentricity e, 
given by Th(e) = r^{e = 0)(1 — 0.16e). We find using 
this formula that the velocity distribution changes only 
by a small amount (less than 2% at 5 km/s) compared 
to when using the Roche lobe radius evaluated at peri- 
center, given by r^(e = 0) above. For the following cal- 
culations, we therefore adopt the Roche lobe evaluated 
at pericenter for an eccentric orbit. 

The radius of each star is found by estimating its ef- 
fective temperature from an isochrone of g iven age t„ in 
the st ellar population synthesis model of iGirardi et al.l 
(2004). This, together with its magnitude, provides an 
estimate of the stellar radius. We denote Mv(m;t g ) 
and R(m; t g ) as the absolute V-band magnitude and 
radius (respectively) of a star of mass m assigned by 
an isochrone of age t g . If the star lies on the horizon- 
tal branch or asymptotic giant branch, instead of using 
its present radius (which may be small) we compare its 
Roche lobe to the largest radius previously attained by 
the star at the end of its red giant phase. 

We shall express our formula in terms of an absolute V- 
band upper magnitude limit Mu m and age t g . Assuming 
the lower magnitude limit to be near the tip of the red 
giant branch, we find the velocity distribution is quite 
insensitive to the exact value of the lower magnitude limit 
because the suppression of binaries due to Roche lobe 
overflow dominates the high-luminosity end of the red 
giant branch. 

Averaging over distributions for all the model param- 
eters and dropping the z subscript for readability, we 
obtain 
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where <d[x] is the Heaviside step function. The variable 
A; is a function of all the other parameters according to 
(j4]), with cf> = <fi(r)). We use a Monte Carlo simulation to 
perform the integration over a grid of log \ v\ values and 
interpolate to find /&(log |t>|). In fig.Q]we plot /h(log \v\) 
for different ages and absolute magnitude limits. This 
figure shows that for velocities > 10 km/s, suppression of 
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Fig. 1. — Distribution of velocities in the center-of-mass frame 
of binary systems, plotted for different absolute magnitude limits 
Mum = 0, 3 and stellar ages t g = 1 Gyr, 10 Gyr. The suppres- 
sion of binaries due to Roche-lobe overflow becomes important for 
velocities > 10 km/s. Except at the turnover point v ?s 10 km/s, 
the distribution behaves locally as a log-normal to good approxi- 
mation. 

binaries due to Rochc-lobc overflow becomes important. 

Asymptotically for large velocities, the velocity distri- 
bution in cq. [6] behaves as a log- normal with dispersion 
similar to that of the period distribution, a\ os p. This can 
be seen as follows: first, at large velocities, the log-normal 
in the integrand is far from its maximum and therefore 
varies slowly in the model parameters (eq. [4j compared 
to their respective probability distributions, provided k 
is not large and negative. Therefore as a first approxima- 
tion we can apply the method of steepest descents and 
find the resulting distribution to be log-normal with a 
mean given by log |u| = (fc — fi\ a gp)- However, the ap- 
proximation of slowly varying log-normal is not strictly 
true, since k becomes large and negative if the mass ra- 
tio q and direction function g e (6,%l>,(p) are close to zero. 
This means that the mean kis in fact a function of v. We 
find, however, that locally k is linear in log \v\ to good 
approximation, with the result that the velocity distri- 
bution still behaves locally as a log-normal but with a 
somewhat different dispersion from o~i og p. We will use 
this to construct a fitting function for the binary velocity 
distribution in section [7] 

3. VELOCITY DISTRIBUTION OF A POPULATION OF 
BINARY STARS 

3.1. Single- and two-epoch likelihood functions 

Suppose that among a population of stars, a fraction 
B of them are in binary systems. Further suppose that 
the velocity distribution for stars not in binary systems 
is Gaussian with dispersion ctq and systemic velocity v. 
The velocity likelihood function will have the following 
form: 



L(v\B,a Q ,v) = (1-B)- 



~(v-v) 2 /2al 



BL b (v\a ,v) 
( ? ) 

where Lf,(v\ao,v) is the likelihood for binary stars. To 
derive the binary likelihood, we note that the component 
of the velocity not due to the binary orbit is v nb = v — v', 



where v' is the line-of-sight component of the velocity in 
the center-of-mass frame of the binary system. To find 
the binary likelihood we therefore average the velocity 
distribution in v n b over the distribution of the binary 
component v' given in eq. [6j 



L b (v\a ,v) = 



2|i/|lnl0 



dv' (8) 



The factor of 2 in the denominator arises from the fact 
that /& (log \ v\) is normalized in log \v\, whereas the like- 
lihood is normalized in v (allowing for positive and neg- 
ative velocities). 

By taking the second moment of the velocity distribu- 
tion in eq. [Tj one obtains the result that 
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(9) 

where a is the measured dispersion and a b is the binary 
dispersion found by taking the second moment of the bi- 
nary velocity distribution fb(v) — /b(log |«|)/2|u| In 10. 
As in the usual case, it can be shown that given a nor- 
mally distributed measurement error with dispersion a m , 
one need only make the replacement Ctq — > Ctq + cr^ in 
the above formulas. 

Next, it is desirable to have a likelihood function for ve- 
locities measured at two different epochs. Since velocity 
changes of order km/s over a timescale of years is entirely 
negligible for nonbinary stars, the most fruitful approach 
is to use a likelihood in the difference Ad between the two 
velocities. Keeping in mind the log-normal behavior of 
the velocity distribution, we write the binary part of the 
likelihood as <7b(log |Au|; At). As with the single-epoch 
velocity distribution /(,(log|u|), we use a Monte Carlo 
simulation to calculate <7f>(log \Av\; At). For each binary 
in the simulation we find Av by evolving the orbital phase 
to its value after a time At and calculating the resulting 
change in velocity. In the absence of measurement error, 
the nonbinaries will have zero change in velocity, so the 
total likelihood can be written as 



L(Av\At, B) = (1 - B)S(Av) 



B 2|AdlnlO (10) 



If there is a normally distributed measurement error, 
the likelihood must be averaged over two Gaussians of 
widths cr m .i and a m ^ for the first and second velocity 
errors, respectively. A little calculation shows this to 
be equivalent to averaging over a single Gaussian with 
dispersion o^e , which is the equivalent 2-epoch measure- 
ment error: 
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Note that the likelihood is identical in form to that 
of eqs. [7]and[8l since in both cases we are averaging the 
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distribution over a Gaussian. Both the single- and multi- 
epoch likelihoods will be put to use in later sections. 

3.2. Threshold fraction of a binary population 

A convenient observable quantity for characterizing a 
binary population is the threshold fraction, defined as the 
fraction F of stars in a sample which exhibit a change 
in radial velocity greater than a threshold Av after a 
time At between measurements. For Aw > 1 km/s, this 
fraction is typically smaller than 0.2, so the threshold 
number (given by n = NF where N is the number of 
stars) follows a Poisson distribution with mean n = NF. 
Therefore the distribution of F is characterized by a sin- 
gle number, the mean threshold fraction F, and the ex- 
pected error can be estimated. In particular, the error in 
F is approximately y/F/N (appendix^. For notational 
simplicity, for the remainder of this paper we will refer 
to the mean threshold fraction F as simply the thresh- 
old fraction F (without the bar), with the understanding 
that the observed threshold fraction will have a Poisson 
scatter about this value. 

Despite its straightforward definition, there are two dif- 
ficulties in measuring the threshold fraction from actual 
data sets. First, often there does not exist a common 
time interval At between measurements in the sample, 
but rather several time intervals for various subsets of 
stars. Furthermore, different velocity measurements have 
their own associated measurement errors and this in turn 
affects the measured value of F. The latter issue can 
be dealt with in an approximate way by using the me- 
dian (or other robust location estimator) of the measure- 
ment error of the sample, in terms of which a2e — 
(cq. Hip . However, both problems can be surmounted 
more rigorously by estimating the error-free threshold 
fraction Fq via a Bayesian or maximum-likelihood ap- 
proach. By using the likelihood in Av defined in eq. 1121 
the threshold fraction at a particular threshold and time 
interval can be estimated even if measurements were 
taken at various epochs — moreover, the inferred thresh- 
old fraction Fq is free of measurement error. This method 
will be demonstrated in sectional 

The threshold fraction without measurement error, 
which we denote by Fq, can be expressed in terms of the 
binary two-epoch velocity distribution <7&(log |Au|; At) 
by taking the integral of eq. QJJ with respect to \Av'\ 
from a threshold Av to oo: 



F {Av\At,B) = B 



Ac 



g6 (log|Ai/|;Af) 
|At/|lnlO 



d\Av'\ (13) 



Note that in the absence of measurement error, the 
threshold fraction Fq scales linearly with the binary frac- 
tion B. The threshold fraction with measurement error 
is likewise obtained by taking the integral of eq. [JS] from 
Av to oo, with the result 



F(Av\At, B, 2e ) = (1 - S)erfc 
'Av - Av' 



Av 



B 



erfc 



V2i 



0~2e 



V2a 2l 
5b (log|Ai/|;At) 



(14) 



2|Au'|lnlO 



d(Av') 



where o~2e is the 2-epoch measurement error given by 



eq. [TTJ Note that in the limit as a2e — > 0, the first term 
goes to zero while the complementary error function in 
the integrand reduces to a step function 20(Ai/ — Av), 
so that eq. [14] reduces to eq. [13] as expected. 

4. CONSTRAINING PROPERTIES OF A BINARY 
POPULATION BY MULTI-EPOCH OBSERVATIONS 

In this section we investigate how properties of a pop- 
ulation of binary stars affect the observed velocity dis- 
tribution measured at two or more epochs. Specifically, 
we explore how our proposed observable, the threshold 
faction F (section |3.2[) . will be affected by changes in the 
underlying binary fraction B, absolute magnitude limit, 
stellar age, size of the measurement error, and time in- 
terval between measurements. We will also demonstrate 
how the binary fraction B can be inferred by a likeli- 
hood analysis, and show how this leads to a better de- 
termination of the threshold fraction F. We first con- 
sider binary models with our fiducial period distribution 
(inferred from the solar neighborhood) and then explore 
how the inferred binary fraction is affected if the assumed 
period distribution parameters are incorrect. Unfortu- 
nately, and as we discuss more fully in the next section, 
the effect of changing the binary fraction B on the ob- 
served binary velocities can be mimicked closely by al- 
tering the assumed distribution of orbital periods (i.e. 
changing the parameters /ii og p, cr\ os p)- While this is 
bad news for any attempt at constraining the underlying 
properties of a galaxy's binary population in full gen- 
erality, it turns out to be good news for correcting the 
observed velocity dispersion for the effects of binary or- 
bital motion, as we will show in section [5] 

First, let us make the rather optimistic assumption 
that the distribution of binary orbital periods is approx- 
imately universal, so that it follows our fiducial choice 
MiogP = 2.23, ciogp = 2.3 (section^. Before launch- 
ing into the full-fledged calculation, one would like to 
estimate how well the fiducial binary fraction B can be 
constrained for a given sample, or conversely, how many 
stars are required to constrain B by a certain amount. To 
simplify matters, let us assume we have a data set where 
the two epochs have the same time interval At between 
them and the same measurement error a m . The equiv- 
alent two-epoch measurement error is then o~2e = cr m v2 
(eq.ril]). 

First, consider the mean threshold fraction of 
the binaries without measurement error, denoted by 
F (Av\At,B = 1) (eq. [13]). In that case the threshold 
fraction scales with the binary fraction B, i.e. is given 
by F (Av\At,B) = BF {Av\At,B = 1). Now consider 
the threshold fraction without binaries, but with a mea- 
surement error a m . The two-epoch measurement error 
is then U2e ~ V^cr m (eq. ITT]) , and the threshold fraction 
F(Av\At, B — 0, <72e) is given by the complementary er- 
ror function (first term in eq. [JJ]with B = 0). In fig. [5] 
we plot the threshold fraction Fo(Av\At, B) produced by 
a Monte Carlo simulation with an absolute magnitude 
limit Mn m = 1 and stellar age t g = 10 Gyr. The thresh- 
old fraction is plotted for different binary fractions, and 
we also plot the threshold fraction from measurement 
error with o~ m — 2 km/s. Near the point of intersec- 
tion Avtaii where F {Av tai i\At, B) = F(Av ta a\At,B = 
0, 02e), the effect of binary stars becomes noticeable over 
the measurement error. Since the Poisson errors are 
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Fig. 2.— Threshold fraction F(Av\ At, B, <r m ), defined as the 
fraction of stars with observed change in velocity greater than 
a threshold A?; after a time interval At between measurements. 
The solid curve has no binaries (B = 0) and a measurement error 
dm = 2 km/s. The other curves are plotted from a Monte Carlo 
simulation for binary fractions B = 1 and 0.5, with At = 1 year 
and no measurement error. The stellar population has an age t g 
= 10 Gyr and the absolute magnitude limit Mi im = 1. Given a 
measurement error <r m , the binary fraction can be constrained for 
thresholds Av > Av tai i, where Av tai i is the point of intersection 
where F(Av\At,B = 0, a m ) = F(Av\At, B, a m = 0). For a given 
binary fraction B, the total threshold fraction without measure- 
ment error is given by B X F(Av\ At, B = 1, a m = 0). 

larger at higher velocity thresholds, to first approxima- 
tion we can say that the error-free threshold fraction Fq 
is best constrained at thresholds near Avtaii- It follows 
that the fiducial binary fraction will be constrained by 
the stars with Av > Avtaii ■ (For a rough approximation, 
one can also use Avt a a ~ 2<72e ~ 2^/2a m .) A little alge- 
bra (see appendix [5]) shows that to constrain the binary 
fraction to within a fractional accuracy of e b , the number 
of stars required is approximately 



N(e b ) 



1 



2B 



F(Av tai i) \ e b 



(15) 



In fig. [3]we graph the approximation formula for differ- 
ent values of B and compare to the 95% confidence inter- 
val in the binary fraction inferred by a Bayesian analysis 
of the simulated data (described later in this section). 
As is evident for the B = 0.7 curve, the approximation 
formula differs for high binary fractions because B > 1 is 
not allowed in the Bayesian analysis. The approximation 
formula is discussed further in appendix [3J_ 

In fig. [4] we plot the threshold fraction F(Av\At, B = 
lj02e = 0) produced by a Monte Carlo simulation with 
an absolute magnitude limit Mu m = 1 and stellar age 
i s = 10 Gyr for different time intervals At. We find that 
for a measurement error a m = 2 km/s, the observable 
threshold fraction steadily increases as At is increased, 
until roughly At — 1 year. This result depends somewhat 
on the mass ratio distribution, since higher mass ratios 
result in higher velocities fo r a given orbital period . If the 
mass ratio distribution in IDuauennov fc Mavorl (|1991[ ) 
is assumed for all periods (as opposed to the uniform 
distribution we adopt for P < 3000 days) , the observable 
threshold fraction increases until roughly At = 2 years. 
In any case, unless the measurement error is smaller than 
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Fig. 3. — Frac tional constraint on the binary fraction B, de- 
fined by eq. 1151 In this plot we a measurement error a m = 2.0 
km/s. For comparison we also plot the 95% confidence interval 
obtained by a Bayesian analysis of simulated two-epoch data from 
a random sample of N stars, averaged over a hundred realizations. 
A uniform prior is assumed for B. 
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Fig. 4. — Threshold fraction F(Av\ At, B, a m ), defined as the 
fraction of stars with observed change in velocity greater than 
a threshold A?; after a time interval At between measurements. 
The solid curve has no binaries (B = 0) and a measurement error 
<Jm = 2 km/s. The other curves are plotted from a Monte Carlo 
simulation for different time intervals At, with no measurement 
error and a binary fraction B = 1. The stellar population has 
an age t g = 10 Gyr and the absolute magnitude limit Mu m = 1. 
Given a measurement error cr m , the binary fraction can be con- 
strained for thresholds Ad > Av ta n, where Avtaii is the point of 
intersection where F(Av\At,B = 0, <x m ) = F(Av\At, B, a m = 0). 
For a given binary fraction B, the total threshold fraction without 
measurement error is given by B X i^Ai^Ai, B = 1, u m = 0). 



2 km/s, little is gained by extending the interval from 1-2 
years to 5 or more years. 

In fig. [5] we plot the threshold fraction for different ab- 
solute magnitude limits Mu m — 0, 3 and stellar ages t g 
= 1 Gyr and 10 Gyr. Extending the magnitude limit to 
fainter magnitudes increases the threshold fraction be- 
cause there is a greater contribution from smaller stars 
with less binary suppression due to Roche-lobe overflow. 
The threshold fraction is also higher for a younger stel- 
lar population because of their larger mass at a given 
stage of stellar evolution, which produces higher orbital 
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Fig. 5. — Threshold fraction F(Av\At,B,a m ), defined as the 
fraction of stars with observed change in velocity greater than a 
threshold An after a time interval At year between measurements. 
Here we have picked At = 1 year, binary fraction B = 1, and no 
measurement error (cr m = 0). We plot the threshold fraction for 
different absolute magnitude limits Mn m = 0, 3 and stellar ages t g 
= 1 Gyr, 10 Gyr. For a given binary fraction B, the total threshold 
fraction is given by F(Av\At, B, a m = 0) = B X F(Av\At,B = 
l,cr m = 0). 
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Fig. 6. — Posterior probability distribution of the binary fraction 
B of a simulated galaxy with binary fraction B = 0.5 and with a 
period distribution characterized by /ii og p = 1.5, <ri og p = 2 (P 
in years). The simulated data sample consisted of 300 stars, each 
with two velocity measurements taken At^ ata = 2 years apart. 
Solid curve is the posterior calculated assuming the fiducial (solar 
neighborhood) model, which is incorrect for this galaxy. Dashed 
curve uses the correct model, with the same period distribution 
parameters /ii og p, o"i og p given above. 

velocities. We find that the threshold fraction changes 
little for ages between 2-10 Gyr; as the stellar age t g 
is reduced from 2 Gyr, however, the threshold fraction 
steadily increases. 

To estimate the binary fraction in a two-epoch sample, 
we use the likelihood function L(Av\At, B, cr 2 e) (eq. IT2)) . 
For the sake of illustration, we analyze a simulated galaxy 
with binary fraction B = 0.5 but with a different period 
distribution from that of the solar neighborhood. We 
choose the period distribution parameters ^i og p = 1.5, 
oiogP = 2 (P in years) for this galaxy. The simulated 
data sample consists of 300 stars, each with two velocity 
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Fig. 7. — Posterior probability distribution of the threshold frac- 
tion F of a simulated galaxy with binary fraction B = 0.5 and with 
a period distribution characterized by fii os p = 1.5, <ri og p = 2 (P 
in years). The threshold fraction F(Av\At) is defined as the frac- 
tion of stars with observed change in velocity greater than a thresh- 
old Av after a time interval At year between measurements. Here 
we chose a threshold Av = 5 km/s and time interval At = 1 year; 
the correct threshold fraction for this galaxy is P ~ 0.05. The sim- 
ulated data sample consisted of 300 stars, each with two velocity 
measurements taken At^ ata = 2 years apart, and a measurement 
error of 2 km/s. Solid curve is the posterior calculated assuming 
the fiducial (solar neighborhood) model, which is incorrect for this 
galaxy. Dashed curve uses the correct model, with the same period 
distribution parameters fi\ og p, oi og p given above. Note that the 
correct threshold fraction can be recovered even if the wrong model 
is assumed (in this case, the fiducial model). 

measurements taken Atdata = 2 years apart, and a mea- 
surement error of 2 km/s. First we assume the fiducial 
model (with n\ os p — 2.23, (J\ gp = 2.3, which is incor- 
rect for this galaxy) and, assuming a uniform prior in 
the binary fraction B, generate a posterior in the binary 
fraction. We then repeat this procedure using the cor- 
rect period distribution parameters /xi og p, cri gp in our 
model, whose values are given above. The resulting pos- 
teriors are plotted in fig. [S] This figure shows that the 
binary fraction B is a highly mo del- dependent quantity, 
and given the unknown nature of the period distribution 
of binaries outside the solar neighborhood, the inferred 
binary fraction must be taken with a grain of salt. How- 
ever, the fiducial binary fraction can still be used as a 
relative indicator of the fraction of observable binaries, 
as long as it is interpreted in reference to the fraction ob- 
served in a binary population following the fiducial (solar 
neighborhood) distributions of orbital parameters. 

Although the inferred binary fraction B is very model- 
dependent, this analysis is still useful in that it leads to a 
better determination of the threshold fraction F. which 
is more directly observable than the binary fraction. To 
see this, we use the Monte Carlo to generate the binary 
threshold fraction of each model, for a threshold Av — 
5 km/s, time interval At — 1 year, and zero measurement 
error, i.e. cT 2e =0 (see eq. [T^|) . We then transform each 
posterior in fig. [B] from B to the threshold fraction F — 
BFf,. The renormalized posteriors P(F) are plotted in 
fig. [71 the correct threshold fraction for this galaxy is 
F 0.05. Note that the correct threshold fraction can 
be recovered even if the wrong model is assumed (in this 
case, the fiducial model). 
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Fig. 8. — Posterior probability distribution for binary fraction of a 
simulated galaxy with 500 stars and binary fraction B = 0.5. Solid 
curve is calculated from three velocity measurements, whereas dot- 
ted curve uses only the first two velocity measurements. v\ and 
V2 were taken 1 year apart, while V2 and 113 were taken 10 years 
apart. 

Since even the stars with velocities smaller than the 
threshold Av are used in the likelihood analysis, the er- 
ror in the threshold fraction F is smaller than if F is 
measured directly, especially for higher velocity thresh- 
olds. The approximate error in the threshold fraction 
F estimated by this technique is derived in appendix [2] 
and given by eq. 1431 Furthermore, the threshold fraction 
at At = 1 year is recovered even though the data was 
taken with a time interval of Atdata = 2 years. More 
generally, the threshold fraction for a specific time in- 
terval can be recovered by the likelihood analysis even 
if the data is taken at various different epochs and with 
different measurement errors. As we will show in sec- 
tion [5j the threshold fraction can be used to correct the 
measured velocity dispersion of a sample for the effect of 
binary motion. 

Finally, one may naturally wonder: how much are 
the constraints improved by including more than two 
epochs in the analysis? To address this question, we 
do a similar calculation on simulated 3-epoch data us- 
ing the Monte Carlo to generate a 3-epoch likelihood, 
L(Av2i, Av 32 |At2i, Ai 32 , B, a^ e ), where the indices refer 
to three velocity measurements v±, V2, and U3, and a^e 
is the 3-epoch measurement error defined analogously to 
eq. QTJ We generate a data sample of 500 stars from a 
simulated galaxy with binary fraction B = 0.5 and the 
period distribution parameters taking their fiducial val- 
ues. For comparison, we generate a posterior P(B) from 
the two-epoch calculation which ignores V3. The results 
are plotted in fig. [8] The velocity measurements v% and 
V2 were taken one year apart, while i>2 and ^3 were taken 
ten years apart. While the most probable inferred binary 
fraction Bfa did improve in this case, the 95% confidence 
limits are only improved by ~ 25%. The fractional im- 
provement in the confidence limits is even less for smaller 
data sets; this is because in a sample of a few hundred 
stars, there is significant scatter in the binary fraction, 
and the inferred binary fraction Bfn has in fact a signif- 
icant probability of becoming worse when more epochs 
are added to the sample. To constrain the binary and/or 



threshold fraction, we therefore find it a more profitable 
strategy to make two-epoch measurements over a larger 
sample of stars, as opposed to adding more repeat mea- 
surements over an existing sample (assuming a similar 
number of overall measurements in either case) . 

A possible complicating factor in the above analysis is 
that selection criteria for making repeat measurements 
can bias the inferred threshold fraction and binary frac- 
tion. If stars whose spectra yield multiple peaks in the 
cross correlation function are singled out for repeat mea- 
surements, the multi-epoch sample may have an inordi- 
nately high binary fraction compared to the overall stel- 
lar population. This selection bias is probably not signif- 
icant in red giant stars due to the typically large differ- 
ence in luminosity between primary and secondary star. 
However, in samples that contain a significant fraction of 
main sequence stars the bias may be more problematic, 
although an upper bound on the binary dispersion can 
still be obtained. 

5. CORRECTING THE OBSERVED VELOCITY 
DISPERSION FROM MULTI-EPOCH DATA 

In the previous section we demonstrated how uncer- 
tainties in the underlying period distribution can ad- 
versely affect our ability to constrain the underlying bi- 
nary fraction from multi-epoch data. Here we demon- 
strate that if our goal is to correct the observed velocity 
dispersion for the effects of binary stars, the degeneracy 
between the period distribution parameters and binary 
fraction is quite useful: regardless of the precise nature of 
the binary population, its effect on the observable thresh- 
old fraction F can be directly related to the associated 
correction in the observed velocity dispersion in a model- 
independent way. 

The important degeneracy arises from the log-normal 
behavior of the binary velocity distribution /(,(log|u|) 
(eq. [6]) . Binary orbital motion along the line of sight of 
order km/s is the most important for the intrinsic disper- 
sions we are interested in. For these velocities, the value 
of log \v\ is far from the mean of the log-normal, which is 
approximately —1 for a magnitude limit Mu m = 1 and 
age t g — 10 Gyr. The exponent of the log-normal is ap- 
proximately linear over the scale of km/s, so we can write 
it as —p — ahx \v\. Therefore the binary part of the veloc- 
ity distribution can be written as f(v) cx Be~@\v\~ 1 ~ a , 
where B is the binary fraction. If the mean binary pe- 
riod /iiogp is varied, the log-normal is offset in log-space 
so to good approximation only /3 changes; therefore the 
velocity distribution f(v) changes by a constant factor 
over the scale of km/s. If the dispersion of the period 
distribution (T\ gp is varied, both the offset /? and the 
slope a change; however, the slope changes by a rela- 
tively small amount for a\ og p ranging from 1-3 (its viable 
range of values; see section [8|), so again the velocity dis- 
tribution changes by an approximately constant factor. 
The important point is if that the parameters a\ og p and 
/iiogp are varied, they change the velocity distribution 
by an amount which is nearly the same over the scale 
of several km/s — in other words, they behave similarly 
as if the binary fraction were changed. This is also true 
of the magnitude limit and stellar age, which effectively 
change the mean of the log-normal and therefore behave 
similarly to n\ og p. We therefore conclude that the pa- 
rameters /iiogp, ciogP, magnitude limit Mu m and stellar 
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age t g are all nearly degenerate with binary fraction over 
the scale of km/s. 

The degeneracy of the period distribution parameters 
with binary fraction also holds for the two-epoch velocity 
distribution (ft, (log | Av\; At), since this also has a log- 
normal form for km/s velocities. By eq.H4| therefore, the 
same degeneracy holds for the threshold fraction F. The 
effect of this degeneracy on the threshold fraction and 
its implications for constraining the binary distribution 
of periods will be explored in further detail in section |8l 

In this section we will consider the threshold fraction 
Fq with a fixed time interval of 1 year and without mea- 
surement error, i.e. Fq = F(Av\At —I year, <72e = 0) 
(section l3.2p . There is no loss of generality in this; as 
we demonstrated in section [4] the threshold fraction for 
any given time interval At can be estimated by a likeli- 
hood analysis even if measurements are taken at various 
different epochs and with various different measurement 
errors. However, if the threshold fraction F is measured 
directly for a fixed time interval, it is necessary to ac- 
count for the effect of measurement error on F; we will 
address this later in the section. 

In the absence of measurement error, by definition 
Fq scales linearly with the binary fraction B (eq. [T3|) . 
Furthermore, because of the near-degeneracy of the pe- 
riod distribution parameters with binary fraction, F also 
scales linearly with /zi og p and a\ og p to good approxima- 
tion over their viable range of values (roughly 1-3 with 
P in years; see section [8|). The essential point is that a 
similar relationship holds for the velocity dispersion if a 
high-velocity cutoff is used, e.g. at v c — 3<r, since the 
degeneracy approximately holds for velocities v < v c . It 
follows that if velocity outlier stars are excluded in de- 
termining velocity dispersion, the extra dispersion due to 
binaries can be determined from the threshold fraction 
Fq with reasonable confidence even if the parameters B, 
WogP, and oiogp are entirely unknown. 

We demonstrate this by simulating galaxies with var- 
ious intrinsic dispersions and characterized by different 
binary populations. The dispersion a is calculated by it- 
eratively discarding stars with velocities larger than 3<r; 
on the first iteration the biweight is used to estimate 
the dispersion, and the dispersion is then calculated on 
every subsequent iteration until all the remaining stars 
have velocities that lie within 3er. In order to make the 
statistical error negligible, we used a very large "sample" 
of 100,000 stars. We also calculate the threshold fraction 
Fo for the same data set, for which we picked a threshold 
Av = 5 km/s and time interval At = 1 year. 

First, we assume the fiducial binary period distribu- 
tion (/iiogp and criogp) and vary the binary fraction from 
B = 0.1 to 1. In fig.[9]we plot the ratio a/ao of measured 
dispersion over the intrinsic dispersion as a function of 
threshold fraction, for galaxies with intrinsic dispersions 
of 4, 7, and 10 km/s. We used an absolute magnitude 
limit Mu m = 3, however the graph remains virtually un- 
changed for other magnitude limits because of the near- 
degeneracy of magnitude limit with binary fraction dis- 
cussed at the beginning of this section. We see that for 
a given intrinsic dispersion, the observed threshold frac- 
tion can be mapped in a one-to-one way to the intrinsic 
dispersion. The relation shown in the graph also holds 
regardless of the age of the stellar population, again be- 
cause of the degeneracy of age with binary fraction. 
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Fig. 9. — Ratio of measured velocity dispersion <r over the intrin- 
sic dispersion cro, plotted with respect to threshold fraction Fq for 
different binary fractions in galaxies of intrinsic dispersions <ro = 5, 
7, and 10 km/s. The measured dispersions were calculated by an in- 
tcrativc 3<r-clipping routine, and the threshold fraction Fq denotes 
the fraction of stars with observed change in velocity greater than 
a threshold Av = 5 km/s after a time At = 1 year between mea- 
surements, assuming zero measurement error. The fiducial period 
distribution (/ii og p = 2.23, cri og p = 2.3, P in years) is assumed, 
and the binary fraction is varied between 0.1 and 1. 

Next we repeat the procedure over a grid of values for 
the parameters B, ^i\ og p, and a\ og p, and for each point 
we plot the ratio a/ao with respect to the threshold frac- 
tion Fo. The results are plotted in fig. [Till again for galax- 
ies with intrinsic dispersions of 4, 7, and 10 km/s. We 
see that for each group, the graph forms a tight relation 
for all but the most extreme values of the period distri- 
bution parameters. In plotting these points we varied B 
from 0.2 to 1, /J,i og p from -1 to 4, and <J\ og p from 0.5 to 
4 (with P in years). The lowermost points of each group 
are the points for with a\ og p has its smallest value of 
0.5, producing only a very small number of short-period 
binaries. The uppermost points are the points for which 
MiogP has its smallest value, so the period distribution is 
shifted toward short periods. For these extreme values, 
the velocity distribution becomes distorted into a dis- 
tinctly non-Gaussian shape so these can be considered 
highly improbable configurations. We have also varied 
the ellipticity distribution parameters e, a e (section [2]) 
and find that the tight correlation in fig.[T0]is unchanged, 
and although the amount of scatter increases slightly, the 
correction still holds to within a few percent accuracy. 

We find that the points plotted in fig.[TU]are well fit by 
a line plus an exponential function, provided that outlier 
points are discarded. To define "outliers" , first we di- 
vide the domain Fo S [0, 0.2] into sections small enough 
so the graph is approximately linear within each section. 
We then further divide each section into two subsections 
and calculate the median and median absolute deviation 
(MAD) of the y-values of the points in each subsection. 
Next we draw lines through the two points defined by 
the median ± twice the MAD of each of the two subsec- 
tions, taking the center of the subsection as their x-value. 
The plotted points that lie outside the region defined by 
these lines represent extreme and highly improbable pe- 
riod distributions, and are therefore discarded. We fit the 
remaining points and repeat the procedure for galaxies 
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* ct = 4 km/s 

* <j Q = 7 km/s 
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Fig. 10. — Ratio of measured velocity dispersion a over the in- 
trinsic dispersion ctq, plotted with respect to threshold fraction Fq 
for different binary populations in galaxies of intrinsic dispersions 
cro = 5, 7, and 10 km/s. The measured dispersions were calculated 
by an interative 3<r-clipping routine, and the threshold fraction 
Fo denotes the fraction of stars with observed change in velocity 
greater than a threshold Av = 5 km/s after a time At = 1 year 
between measurements, assuming zero measurement error. Each 
point represents a different binary population with its own binary 
fraction and period distribution; we plotted the points over a grid 
of values, with binary fraction B ranging from 0.2 to 1, niogP 
from -1 to 4 (in log(P/year)) and criogp from 0.5 to 4. We show 
at the end of section [5] how Fo can be inferred from observations 
accounting for measurement errors. 

of dispersions ranging from 3-12 km/s. 

The plots in fig. [TU] are not directly applicable to real 
data because each graph was plotted for a fixed intrinsic 
dispersion Co, which is unknown (and is in fact what we 
are attempting to calculate!). We therefore use our fits 
together with a root-finding procedure, interpolating the 
fitting parameters in cto, to draw similar graphs at fixed 
values of a. A few resulting curves are plotted in fig. 1111 
Again, we find these curves are well fit by a line plus 
exponential, 



(TO 



o(cr) + b(a)F + c(a) 



,F /0.1 



1 



(16) 



where Fo is the threshold fraction at 5 km/s. We also 
find fitting functions for the parameters a(cr), b(a), and 
c(ct) which fit well for dispersions ct ranging from 4 km/s 
to 10 km/s. Defining Act = a — 4 km/s, we find: 

a(CT)=0.988e-°- 0007Aff (17) 
b(a) = 0.576 - 0.08Act + 0.772 (l - e" 01Acr ) (18) 
c(CT)=0.043e-°- 247ACT (19) 

These formulas hold for any magnitude limit and stel- 
lar age, and the threshold fraction Fq refers here to a 
velocity threshold of 5 km/s, time interval of 1 year, and 
zero measurement error. 

How are these formulas adjusted if a different veloc- 
ity threshold is desired? Ideally, one should use the 
smallest possible threshold that is not significantly af- 
fected by measurement error-this will include the most 
stars and therefore have a smaller scatter compared to 
higher thresholds. To use a different threshold Av, 
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Fig. 11. — Ratio of measured velocity dispersion a over the in- 
trinsic dispersion <ro, plotted with respect to threshold fraction Fq 
for fixed measured dispersions a. Thes e curves were found by fit- 
ting graphs like those shown in fig. 1101 then transforming to fixed 
val ues of the measured dispersion <r. The measured dispersions in 
fig- HOI were calculated by an interative 3<r-clipping routine, and 
the threshold fraction Fo denotes the fraction of stars with change 
in velocity greater than a threshold Av = 5 km/s after a time 
At = 1 year between measurements, with zero measurement er- 
ror. We show at the end of section l5l how Fo can be inferred from 
observations accounting for measurement errors. 

the value of Fq in eq. [16] must be scaled by the ratio 
Fq(Av) / Fq(5 km/s). This ratio can be calculated by us- 
ing Monte Carlo realizations to plot the threshold frac- 
tion as a function of threshold, which we will do in section 
[4] (figs. [2] [4j. Again, the degeneracy of magnitude limit 
and period distribution parameters with binary fraction 
ensures that this ratio is virtually independent of the 
model parameters and magnitude limit, provided one 
does not transform to thresholds that are too high (> 
10 km/s). We find that for thresholds in the range 4 
km/s < Av < 10 km/s, the ratio F (Av)/F (5 km/s) 
can be fit by the function 



Fg(Av) 

F (5 km/s) 



be -±v/Av s 



(20) 



where the best-fit parameters are a — 0.0725, b — 1.897, 
and Av s = 6.947 km/s. Thus to find the correction to 
the dispersion in terms of a given velocity threshold At), 
one substitutes eq. [20] into eq. [16] so that the fit is in 
terms of F Q (Av). 

The sampling error in the intrinsic dispersion cto deter- 
mined by this procedure can be estimated by noting that 
for a two-epoch sample of N^e stars, the Poisson error in 
the threshold fraction is 6F fa y / F n (Av)/N2 e - Let us as- 
sume the dispersion is measured in a larger single-epoch 
sample of N stars; then we can make the approximation 
that the sampling errors in a and F are weakly correlated 
so they add in quadrature. Propagating the error using 
eq. U6l gives the sampling error in the intrinsic dispersion, 



Sao 



2N 



a d(a/a ) 



ct dF n 



F (Av) 



N, 



2e 



(21) 



where the second term in eq. [21] is the two-epoch sam- 
pling error. For example, if the single-epoch sample con- 
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tains N = 500 stars and the measured threshold fraction 
is Fq = 0.1 in a two-epoch subset of A^e = 100 stars, 
the intrinsic dispersion uq can be determined to within 
« 5% for dispersions a > 4 km/s. Eq. |2"T"1 was tested 
with Monte Carlo simulations and found to be accurate 
to within sa 0.8% for N = A^e! in the above example, 
the formula is accurate to within w 0.5%. However, the 
fractional error in <jq cannot be made smaller than 
0.5% because of the inherent uncertainty in the binary 
population model represented by the width of the plots 
in fig. [lOl We find that for an observed threshold frac- 
tion F — 0.1, the two-epoch sampling error is smaller 
than the single-epoch error unless N/N2 e > 5. Given a 
measured threshold fraction F(Av), the two-epoch sam- 
pling error is larger relative to the single-epoch error for 
smaller measured dispersions. 

Finally, we address the issue of measurement error. 
As mentioned above, the best approach is to estimate 
the error-free threshold fraction Fq by a Bayesian or 
maximum- likelihood approach (see section [4]), in which 
case measurement error need not be considered here. 
However, if F is calculated directly from the data, then 
measurement error must be considered in the above for- 
mulas. We find that given a 2-epoch measurement er- 
ror cr 2e (eq. QT]), the threshold fraction F(Av\<j2e) is re- 
lated to the threshold fraction without measurement er- 
ror Fq(Av) by a linear transformation — yet another con- 
sequence of the degeneracy between binary fraction and 
the period distribution parameters (see appendix [1] for a 
derivation of this result). The transformation takes the 
form: 



F(Av\a 2e ) ~ erfc 



Av 



&2e 



Av 

<?2e 



Fo(Av) (22) 



This result (together with eq. 120]) has been tested by 
using the Monte Carlo simulation to produce plots sim- 
ilar to fig. [10] for different velocity thresholds and mea- 
surement errors. The approximate analytic form of 
can also be derived (see appendix [TJ); using this together 
with the Monte Carlo plots to map f3, we find that j3 is 
well fit by the following function (appendix [TJ : 



(3{x) = (1 + ae 



1 — k ■ erfc 



V2 



(23) 



where x = Av/o^e and the best-fit parameters are a = 
0.05, x s = 5, and R = 1.3. Thus to find the correction 
to the dispersion with a given measurement error, one 
substitutes the error-free threshold fraction Fq in terms 
of F(Av\a 2e ) (given by eq. [22J into eq. Q£l The effect 
of the measurement error on the dispersion must also 
be taken into account by making the substitution a 2 — 
a meas ~ °m ec l- DJH where o m eas is the measured 
dispersion and a m is the measurement error. 

In summary, the velocity dispersion of a dwarf 
spheroidal sample can be corrected for binaries by the 
following method: 

1. Measure the threshold fraction F for a particu- 
lar threshold velocity and a time interval At = 1 year. 
This can be done in two ways: the threshold fraction can 
be measured directly, in which case one should use the 
smallest possible threshold that is not unduly affected 



by measurement error; this is approximately Av « 2a 2e 
where <72 e = v2<7 m is the median two-epoch measure- 
ment error. An alternative (and more rigorous) ap- 
proach is to estimate the measurement error-free thresh- 
old fraction Fq by a likelihood analysis. This procedure 
is demonstrated in section [4j 

2. Measure the velocity dispersion a meas of the sample 
by an iterative 3er-clipping routine. Correct the disper- 
sion for measurement error to find the error-free disper- 
sion a. 

3. If the chosen velocity threshold is different from 5 
km/s, scale the threshold fraction Fo(5 km/s) in eq. IT"6l 
according to eq. [2TJ1 

4. If the threshold fraction is measured directly, one 
must use eq.[22]to relate the error- free threshold fraction 
Fq(Av) in eq.[16]to the threshold fraction with measure- 
ment error, F{Av\a2e)- 

5. After substituting Fq(5 km/s) in terms of F(Av\a2e) 
(given by steps 3 and 4 above) into eq. [121 use eq. [TH] 
together with eqs. lTTlfTIjl to find the intrinsic velocity dis- 
persion (To. The sampling error in <7q can be estimated 
from eq. [21] 

Finally, from fig. [10] we can estimate an approximate 
upper bound to the dispersion introduced by binaries in 
dwarf spheroidals. First we note that given an abso- 
lute magnitude limit Mu m < 1 and a stellar population 
older than 1 Gyr, and if the distribution of binary or- 
bital parameters mirrors that of the solar neighborhood, 
the threshold fraction Fq cannot be larger than 0.12 even 
if the binary fraction is 1 (fig. [5]). By analyzing multi- 
epoch data in the Fornax, Carina, Scul ptor, and Sextans 
dwarf spheroidals ijWalker et al.l I2OO90 we find they all 
have threshold fractions smaller than 0.15, and only For- 
nax has F > 0.1. Assuming this is generally the case 
even for ultra-faint dwarf spheroidals with intrinsic dis- 
persions greater than 4 km/s, we conclude from fig. [10] 
that the measured velocity dispersions of these galaxies 
are unlikely to be inflated by more than 20%. However, 
even if the dispersion of a particular galaxy is inflated 
by more than 20%, the correction due to binaries can be 
readily discerned by making repeat measurements and 
applying the method outlined above. 

6. BAYESIAN ANALYSIS OF SINGLE-EPOCH VELOCITY 

DATA 

In this section we discuss the problem of constrain- 
ing properties of a binary population with radial veloc- 
ity measurements taken at a single epoch. We approach 
this by fitting the likelihood for binary stars (eq. [7J to 
a simulated data sample generated from a Monte Carlo 
simulation for galaxies with dispersions 4 km/s and 10 
km/s. Both galaxies were given a binary fraction B = 0.5 
and period distribution parameters equal to their fidu- 
cial values /Ltiogp = 2.23, <J\ og p = 2.3 (P in years). The 
simulated velocities were generated with a measurement 
error of 2 km/s. Using eq. [7] together with Bayes' The- 
orem, marginal posterior distributions in the dispersion 
a and binary fraction B can be obtained. The maxima 
of these distributions are taken as the best-fit values a fa 
and Bf it . 

In the absence of any constraint on the binary fraction 
from multi-epoch data, we assume a uniform prior in the 
binary fraction. However, if repeat measurements do ex- 
ist for some subset of the data sample, these can be used 
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binary fraction (B) binary fraction (B) 

(a) (T rea i = 4 km/s, 100 stars (b) cr rea i = 4 km/s, 500 stars 

Fig. 12. — Best-fit values for the dispersion afn in a simulated galaxy with dispersion <r rea i = 4 km/s and different binary fractions B. 
The solid vertical line gives the error bar, equal to one standard deviation in <jf it values calculated in 300 random realizations, which is 
similar for all points on the graph. Solid line uses a uniform prior in B, while dashed line uses a prior P(B) obtained from multi-epoch 
observations of the stars. Dot-dashed line is calculated from the biweight robust estimator. We adopted a measurement error of 2 km/s 
for all stars. 
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(a) a rea i = 10 km/s, 100 stars (b) oy ea ; = 10 km/s, 500 stars 



Fig. 13. — Best-fit values for the dispersion Of it in a simulated galaxy with dispersion cr rea i = 10 km/s and different binary fractions B. 
Error bar gives one standard deviation in crj it values calculated in 300 random realizations, which is similar for all points on the graph. 
Solid line uses a uniform prior in B, while dashed line uses a prior P(B) obtained from multi-epoch observations of the stars. Dot-dashed 
line is calculated from the biweight robust estimator. We adopted a measurement error of 2 km/s for all stars. 



to better constrain the binary fraction. The most rigor- 
ous approach would use a joint likelihood Av\At), 
which can be generated from the Monte Carlo. This 
would have the advantage that individual binary stars 
which are inflating the dispersion and exhibit a large 
velocity change Av would be recognized as such, and 
weighted accordingly. While this method may be neces- 
sary for obtaining constraints in small data sets of less 
than 100 stars, for larger samples we can adopt a sim- 
pler approach: first derive a posterior P{B) in binary 
fraction by the multi-epoch analysis outlined in section 
IH and subsequently take P(B) as a prior in B for the 
single-epoch analysis. The usefulness of this method de- 
pends critically on the size of the multi-epoch sample, 
as this determines the constraint on binary fraction (see 
fig.©. 

To evaluate this method, it is important to know how 
the best-fit dispersion a fa obtained by this procedure 
may differ from the intrinsic dispersion a rea i due to small 



number statistics. To investigate this, we analyze sam- 
ples consisting of 100 and 500 stars and repeat the pro- 
cedure over 300 randomly generated realizations of each 
sample. The range of <jfu values obtained for galaxies 
with dispersions of 4 km/s and 10 km/s are plotted in 
fig- H21 and fig. [HI respectively. Also included are the 
results if the prior in B is taken from multi-epoch obser- 
vations of the same stars in the sample, in the manner 
outlined above. For comparison, we also plot the disper- 
sion obtained from the biweight robust estimator, which 
is roughly comparable to the dispersion obtained by us- 
ing a using a 3<7-clipping routine. This plot shows the bi- 
weight is biased to large dispersions by an amount which 
depends on the binary fraction; furthermore, the bias for 
a 4 km/s dispersion is much larger compared to the er- 
ror bars than for a 10 km/s dispersion. By comparison, 
the best- fit dispersion a ja using a uniform prior in B is 
biased by a smaller amount which depends on the true 
binary fraction, up to 10%. However, if the prior P{B) 
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is calculated from multi-epoch samples of as few as 100 
stars, the bias is almost entirely removed. Thus, even a 
fairly weak constraint on the binary fraction is sufficient 
to remove most of the bias. 

Unfortunately, the likelihood analysis outlined above 
is of limited utility in actual data sets because of con- 
tamination by foreground Milky Way stars. The usual 
criteria to determine membership of a star include its lo- 
cation on the color-magnitude diagram, metallicity, and 
radial velocity. Outliers in the distribution of these vari- 
ables are either excluded, or else weighted by a low mem- 
bership pr obability assigned according to a specific algo- 
rithm (cf. IWalker et al.ll200"9l ). However, among the stars 
with single-epoch measurements, only those stars which 
lie on the high-velocity tail of the velocity distribution 
can help constrain the binary fraction. If even a small 
number of high-velocity binary stars are excluded from 
the data or weighted by a low membership probability, 
the inferred binary fraction will be biased to low values. 

In principle this problem can be resolved by including 
the velocity distribution of foreground Milky Way stars 
in the likelihood analysis, obtained b y a kinematic mode l 
of the disk and bulge components (|Robin et al.l 12003). 
The success of this method will depend critically on the 
degree of overlap between the two distributions, which is 
partly determined by the systemic velocity of the back- 
ground galaxy. It is also essential that binarity in the 
Milky Way is accounted for, since this adds a substan- 
tial high-velocity tail to the foreground velocity distribu- 
tion. Accounting for binarity in the likelihoods of both 
the foreground and background stars may also lead to 
improved membership probabilities when c ombined with 
the ex pectation maximation algorithm of Walk er~et al.l 
(2009). Including both binarity and foreground Milky 
Way stars in a likelihood analysis is certainly of in- 
terest for obtaining the best possible constraints from 
single-epoch velocity data, but is beyond the scope of 
the present work. 

7. FITTING FUNCTION FOR BINARY LIKELIHOOD 

For the purpose of making analytic calculations as well 
as doing a likelihood analysis, it is useful to have a fit- 
ting function for the binary velocity distribution. This 
can be found by investigating the exponent of the binary 
velocity distribution. We find that the distribution in- 
deed follows a log-normal but the dispersion changes at 
a velocity scale vq due to the suppression of binaries by 
Roche lobe overflow. This suggests the following fitting 
function for the binary distribution of velocities: 



/ b (io g M) = (l-AOvQogM) 
-l 



Af exp ■ 



2 <P 



[Wo g p+3G(logM)n, (24) 



As we discuss below, all the fitting parameters are con- 
tained in the function Q (log \ v\): 



s(iogM) 



a + s log \v\, 
a + s' log 



M < V 
+ slogvo, M > V . 



(25) 

For velocities \v\ < vq, the log- normal dispersion is 



given by cr log | „ | = F . The destruction of binaries due 
to Roche-lobe overflow becomes important for velocities 
|u| > Vq, for which the log- normal dispersion steepens to 

_ CT log P 

°logM — 3s' ■ 

The normalizing factor Af (given by eq. [57] below) is 
not a fitting parameter, but rather is determined by the 
other fitting parameters along with /ii og p and ai og p . The 
function T>, whose parameters are fixed, is given by 



V{\og\v\ 



-(log\v\-ef/2a 2 , 



(26) 



and is effectively a smoothed S- function at very small v. 
With v in units of AU/year, we find e = —3 and as = 1 
works well with only slight deformation of the distribu- 
tion at small v, however the exact values of these param- 
eters are unimportant so long as the width and mean of 
T> is sufficiently small. Although stars at such small ve- 
locities give a negligible contribution to the dispersion, 
including the smoothed (5-function is crucial as it gives 
the appropriate weight to the small-velocity stars. 
The normalizing factor Af in eq. [24] is given by 



27rfj Lp 



»7(MiogP,0-io g p)' 



(27) 



1 l 

??O«iogP,0-iogp) = ^-erfc[-a;(s)] + — erfc[x(s')], (28) 



x{s) 



3s log \vq \ + 3a + ^iogp 

V^OlogP 



(29) 
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0.210 


0.888 


1.920 


2.445 


0.878 


3.460 


1 


0.240 


0.938 


1.770 


2.389 


0.836 


3.448 


2 


0.195 


0.835 


1.467 


2.430 


0.937 


4.897 


3 


0.152 


0.737 


1.299 


2.470 


1.056 


6.986 



TABLE 1 

Fitting parameters a, s, s' , vo for binary velocity 

DISTRIBUTION (EQ. l2"4l FOR DIFFERENT ABSOLUTE MAGNITUDE 
LIMITS My . WE ALSO TABULATE THE NORMALIZATION JV AND 

BINARY - DISPERSION a b ASSUMING THE FIDUCIAL PERIOD 
DISTRIBUTION PARAMETERS £tl og P = 2.23, (7 log p = 2.3 (P IN 
YEARS). 



The fitting parameters are tabulated in table ?? for 
different absolute magnitude limits Mn m ; for magnitude 
limits in between the tabulated values, interpolation may 
be used to find the fitting parameters. We also tabulate 
the normalization Af and binary dispersion at assuming 
the fiducial period distribution parameters ni og p = 2.23, 
ciogP = 2.3 (P in years). While these parameter values 
were chosen with v in units of AU/year, converting to 
other units simply translates /f,(log \v\) along the log \v\- 
axis by an amount A log \ v\ — logK, where k is the con- 
version factor from AU/year. The conversion factor to 
km/s is k — 4.741. 

The function /f,(log|t>|) in eq. [M]can be conveniently 
used in place of eq. [6] which is computationally difficult. 
Likewise, it can be substituted in eq. [8] and integrated 
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numerically to calculate the likelihood of binary stars in 
a galaxy of given intrinsic dispersion a. 

If the period distribution parameters /Ui og p and a\ og p 
are varied from their fiducial values, the fitting parame- 
ters remain unchanged to first order (although the nor- 
malizing factor changes according to eqs. I27l- l2"9"|) . Thus 
we have the additional benefit that the binary dispersion 
can be calculated and its dependence on the period dis- 
tribution parameters is seen explicitly. Analytic formulas 
for the binary dispersion are given in appendix [3l Here 
we shall simplify matters by taking a log-normal distri- 
bution so that s = s'; in this case, we find the approxi- 
mation dispersion is proportional to 



a b oc exp 



/ oiogp In 10 



3s 



exp 



, (WogP + 3a) In 10 



3.s 



( 3 °) 

Thus the variation of Ub with the period parameters 
is quite distinct from its variation with B, which is lin- 
ear (eq. I!]). If one finds the dispersion up to a cutoff 
velocity (of order AU/year), however, extra factors are 
introduced (see eq. [50] in appendix [3]) which mitigate the 
exponential dependence with <7i og p and ^i og p and nearly 
linearize them over their viable range of values (roughly 
from 0-4 with P in years; see section [8]). Thus if one 
chooses a cutoff v c — ca where a is the measured disper- 
sion (e.g. c = 3), the dispersion is approximately linear 
in the period parameters, consistent with the expected 
near-degeneracy with binary fraction. This reinforces the 
very useful observation that the effect of binaries on the 
measured dispersion of a galaxy can be discerned from 
multi-epoch data in a model-independent way. 

8. CONSTRAINING THE DISTRIBUTION OF PERIODS 

Here we address the question of what form the dis- 
tribution of periods might take in regions beyond the 
solar neighborhood, and whether this can be constrained 
by radial velocity data. Simulations of star f ormation 
(jMachida et alJl2009L lTohlmell2002L lBatell2000( ) indicate 
that the statistical properties of binary systems are de- 
termined during star formation via turbulent fragmenta- 
tion of a rotating gas cloud. The distributions in the 
orbital parameters generally undergo little subsequent 
modification by collisional processes, except in the high- 
density reg ions found at the centers of globular clusters 
(|Hut et al.111991 iPrvor et al.lll988| ). This suggests that 
the distribution of periods in dwarf sphcroidals and dwarf 
irregular galaxies may be of a similar form to that found 
in the solar neighborhood. We shall assume the log- 
normal is a fair approximation to the period distribu- 
tion; it is reasonable to ask, however, to what extent its 
mean /ii og p and dispersion a\ og p may be expected to dif- 
fer from that of the solar neighborhood. Because of the 
difficulty of simulating binary star formation, at present 
we have an incomplete picture of how these parameters 
might vary depending on the star formation history of 
each galaxy. However, the following points can be made. 

A semi- empiric a l mod el of isolated binary star for- 
mation by iFisherl (|2004f ) yielded values of a\ og p within 
the range 1.6-2.1, depending on the star formation ef- 
ficiency of the initial gas cloud. Observations of pre- 
main sequence stars in Milky Way stellar associations 
also show more peaked distributions than a\ og p = 2.3 
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Fig. 14. — The behavior of the threshold fraction F(Av) if the 
width of the period distribution criogp (P in years) is varied, plot- 
ted for different velocity thresholds Av. The y-axis is given by 
/ = F/Fo, where _Fo is the fiducial threshold fraction assum- 
ing the value Ci og p = 2.3 observed in the solar neighborhood 
HDuquennov fc Mayor 199j). 
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Fig. 15. — The behavior of the threshold fraction F(Av) if the 
mean of the period distribution /ii og p (P in years) is varied, plot- 
ted for different velocity thresholds Av. The y-axis is given by 
/ = F/Fo, where _Fo is the fiducial threshold fraction assum- 
ing the value rtogP = 2.23 observed in the solar neighborhood 
HDuquennov fc Mayor 199j). 



(|Brandner fc Koehlerl I1998D . This suggests that the 
wider distribution observed in the solar neighborhood 
may be formed by a superposition of more sharply peaked 
binary distributions resulting from various star-forming 
environments. Accordingly, to be conservative in this pa- 
per we have considered values of ai og p € (0.5,4). Like- 
wise we have considered values of m gP € (—1,4) (P 
in years). The low-^i og p, high-cri og p boundaries of these 
intervals describe velocity distributions with a highly dis- 
torted, non-Gaussian shape and therefore can be consid- 
ered unlikely. 

Supposing that u\ n „ p and cri og p may vary by the 
amount suggested bv lFisherl (|2004D . can these param- 
eters be estimated independently of binary fraction? To 
answer this question, we reconsider the threshold frac- 
tion F{l\v\/At), defined as the fraction of stars with 
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an observed change in velocity greater than a thresh- 
old Av after a time At between measurements. For 
each threshold, let us define the ratio / = F/Fq where 
Fq is the fiducial threshold fraction obtained by setting 
o'logP = 2.3, fJ>iogP = 2.23. Given a time interval of 1 
year and picking several thresholds Aw, we plot / as a 
function of ai og p in fig. [TH and /ii og p in fig. [15] The simi- 
larity in slope among the different curves give a measure 
of degeneracy; if the curves were identical, the param- 
eters would be completely degenerate with the binary 
fraction. However, if the difference between the curves 
is smaller than the statistical error bars, the parameters 
are effectively degenerate with B and we cannot hope 
to distinguish between them. The error bars on / at a 
particular Av are given by (compare eq. [43] in appendix 
12): 



2/ 



y/NF(Av tai i\At) 



(31) 



where e/ is the 95% confidence limit in /, N is the num- 
ber of stars and Avt a a is defined as before (for a rough 
approximation, given a measurement error <r m one can 

take Avtaii ~ 2c2e where <72e ~ V^ffm- For At = 1 
year and a measurement error of 2 km/s, fig. Ogives 
F(Av t aii) « 0.1 x B. If B « 0.5, this gives e f ps 9f/y/N. 
Assuming a sample of N = 100 stars, we have £/ ~ 0.9, 
far too large to distinguish between the different curves 
in figs. [T4l and [T5l Thus to within the statistical error in 
a sample of a few hundred stars, a\ og p and /ii og p are de- 
generate with the binary fraction over nearly the entire 
range of figs. [T4l and ITSl only with N > 1000 are the error 
bars small enough to break this degeneracy. We conclude 
that multi-epoch samples of ps 1000 stars or more are re- 
quired in order to constrain the period distribution of a 
population of binary systems. 

9. CONCLUSION 

We have demonstrated a procedure to estimate the in- 
trinsic velocity dispersion of dwarf spheroidal galaxies 
more precisely than in previous studies by accounting 
for the effect of binary orbital motion. By measuring 
the threshold fraction of a stellar sample (defined in sec- 
tion 13. 2p , the correction to the velocity dispersion due 
to binary motion can be estimated; fitting functions are 
provided for this method (eqs. [T6l fT7lfT9|) . We have also 
demonstrated a method to estimate the threshold frac- 
tion, together with other properties of the binary pop- 
ulation, more rigorously from multi-epoch data via a 
Bayesian or maximum likelihood approach. We conclude 
with the following points: 



1. The dispersions of dwarf spheroidal galaxies be- 
tween 4-10 km/s are unlikely to be inflated by more than 
20% by binaries (fig. [TU]). This conclusion can be made 
with confidence because — as we showed in section[5] — the 
correction to the dispersion holds independently of the 
model used to describe the binary population, provided 
that the dispersion is calculated using a high-velocity 
cutoff (e.g., by a 3<7-clipping routine). This is fortunate 
since, as we showed in section [5] the binary fraction and 
other properties of the binary population (e.g., period 
distribution) are very difficult to constrain independently 
of each other with samples of only a few hundred stars. 

2. With a measurement error on the order of 1 km/s, 
we find that 1-2 years is an optimal interval between mea- 
surements for measuring the threshold fraction, since the 
fraction of stars with a measureable change in velocity 
does not rise significantly when the interval is extended 
beyond 2 years (fig. [4]). Furthermore, to constrain the 
binary fraction and other properties of the binary popu- 
lation, it is a more profitable strategy to make two-epoch 
measurements over a larger sample of stars, as opposed 
to adding more repeat measurements to an already exist- 
ing two-epoch sample. This is necessary to overcome the 
large scatter in the binary fraction of samples with less 
than a few hundred stars. We also find that multi-epoch 
samples of more than 1000 stars would be required in 
order to constrain the distribution of periods of a binary 
population independently of the binary fraction. 

3. While the binary population can in principle be 
constrained by single-epoch data, in practice this is very 
difficult because of contamination by non-member stars. 
Since an outlier in the velocity distribution cannot be 
verified as a binary star without multiple epoch mea- 
surements, it may be erroneously labeled a non-member 
star and excluded from the data sample (or weighted 
by a low membership probability). This would bias the 
estimated binary fraction to low values, resulting in an 
inflated dispersion estimate. Even in single-epoch sam- 
ples, however, this problem may be overcome by includ- 
ing a likelihood for the foreground Milky Way stars in a 
Bayesian analysis. This analysis can also be combined 
with multi-epoch data to provide better constraints, and 
in principle would extract the best constraints from both 
single- and multi-epoch velocity data. 
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1. APPENDIX: CALCULATING THE THRESHOLD FRACTION FOR DIFFERENT MEASUREMENT ERRORS 

In section[5]we showed how to correct the dispersion (eq.lTB]) by using the threshold fraction F, defined as the fraction 
of stars with an observed change in velocity greater than a threshold Av after a time At between measurements. If the 
threshold fraction is calculated directly from the data, then measurement error must be considered in eq. [16] While 
the measured dispersion can be easily corrected for measurement error according to cr^eus — a 2 — er^ (where a meas is 
the measured dispersion and a m is the measurement error), correcting the threshold fraction is less trivial. To correct 
the threshold fraction for measurement error, first we note that the degeneracy of the period parameters with binary 
fraction means that the velocity distribution can be approximately split into two parts: a small-w part which acts 
effectively like a 8- function similar to that in eq. 1101 and a large- v tail which scales linearly with B, fJ,\ og p, and <7i og p. 
Exactly where to do the "splitting" is dictated largely by the measurement error, as the tail becomes prominent only 
at velocities beyond Av w a 2e . We therefore split the function at Avt a a — l< J 2e where 7 is a proportionality constant 
with a very weak (if any) dependence on a 2e . The threshold fraction can then be written approximately as follows: 



where 



F(Av\At, B, a 2e ) « (1 - A/>rfc 
'Av - Av' 



Av 



J{Av\a 2e 



erfc 



9b (log I Av' I) 



lAu'llnlO 



BJ(Av\a 2e ) 
d(Av'), 



(32) 
(33) 



Af = J BJ( 7 a 2e |c7 2e ). (34) 

The normalization factor Af here has replaced B in eq. 1 141 and varies linearly with B, /ii og p and u\ 0& p. Now as long 
as Av is approximately equal to or larger than <72 e , the error function in the integrand of eq. [33J is approximately a 
smoothed step function. For Av / a 2e ^ 1 we can therefore make the approximation 

J(Av\a 2e ) « a (^j F (Av) (35) 

where Fq(Av) = F(Av\a 2e = 0) denotes the threshold fraction without measurement error. Note that in the limit 
f2e (Av/a 2e — > 00), the integrand of eq. 1331 becomes a step function so that a — > 1. On the other end, as Av/a 2e 
is made smaller, the integral in eq. [33J includes more of the central peak so a becomes larger. Using eq. 1351 we can 
write Af = kFq(Av) where 

k = " (7)J ? (7g2e) (36) 



Since Fq is the tail of a log-normal distribution, over the scale of several km/s we have Fq(Av) ps (Aw) ™ where n 
is an exponent close to 1. Therefore k can be written as a function of Av/a 2e : 



(a) 



,C2, J ,, [ All 1 • 17 2, 

Substituting this result into eq. [32] we have 
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F(Av\a 2e ) ~ erfc 

where 



At 



v/2, 



0"2e 



C2e 



F (A«) (38) 



Eq. [38] is the key result: it means that the threshold fraction i* 1 with a given measurement error is related to the 
measurement error-free value Fq by a linear transformation, and the only extra information required to make this 

transformation is the function j3 (^f") • As a check, taking the limit a 2e ->0we have a — > 1 and k — > 0, so that /3 — > 1 

as it should. We use the Monte Carlo simulation to plot f3 as a function of Av/a 2e for different velocity thresholds, 
and find the function (3 is nearly the same regardless of threshold, which justifies the approximations taken to reach 
eq. We find that for Av / a 2e ^ 1 , the parameter k shows only slight variation over the range for which the error 
function is non-negligible. We also find that a is well approximated by an exponential, so that /? is well fit by the 
following function: 

P(x) = (l + ae~*) jl-K-erfcf-^=H (40) 

where x = Av/a 2e and the best-fit parameters are a = 0.05, x s = 5, and K = 1.3. 

2. APPENDIX: FORMULA FOR NUMBER OF STARS REQUIRED TO CONSTRAIN THE BINARY FRACTION 

For a set of N stellar velocities (with TV > 100) measured at two different epochs separated by time At, consider 
the number of stars n with change in velocity greater than some threshold value Av. If n were to be measured over 
many random realizations drawn from a particular galaxy, it would follow a Poisson distribution. For the time being, 
let us take the measurement error a 2e to bezero. The mean number h is related to the mean threshold fraction 
of binaries Fb(Av\ At, a 2e — 0) by n = NBFb where N is the total number of stars. If we pick Av small enough 
such that this number is larger than 10, the Poisson distribution is approximately Gaussian with standard deviation 

yj NBF b . Since the fraction of stars is F — n/N, we obtain the error in the measured threshold fraction F: 



j F(Av\At,B,a 2e =0) 
= V n (41) 

where we also used the relation F = BF b . Now if a 2e is nonzero, the smallest value of Av which is largely unaffected 
by the measurement error will occur roughly at the point Av ta u where F(Av ta u \ At, B, a 2e = 0) = F(Av ta u \ At, B = 
0,a 2e ); this is where the "tail" in the distribution due to binaries begins to dominate. Therefore we pick this point 
as giving the best constraint on B. (For a rough approximation, one can also use Avt a u ~ 2&2e-) If we measure the 
fraction of stars with change in velocity greater than Avt a u, given by F = n/N, our "best-fit" binary fraction b is then 
defined by F — bFi,(Avtaii\At, B, a 2e — 0). Combining this with F = BF^ and eq. HH gives the standard deviation of 
the best-fit binary fraction: 



° b \j N x F b (Av tail \At,a 2e = 0) (42) 

Now calling the 95% confidence limit e b = 2a b and solving for N yields eq. Qj)l The approximation becomes less 
accurate as B tends towards very small (close to zero) and large (close to one) values. If B is sufficiently close to zero 
such that n is less than 10, the Gaussian limit no longer holds; in that limit, the uncertainty will be larger than that 
given in eq. HU If B is close to 1, the width of the Gaussian is larger than the true uncertainty since i3-values greater 
than 1 are not allowed. In fig. [3] the approximation formula is graphed as a function of N for several values of B and 
compared to the 95% confidence interval of the posterior P{B) obtained from a Bayesian analysis of the simulated 
data, averaged over a hundred realizations. 

Using the best-fit binary fraction obtained by the threshold fraction at Avt a a , we can find_the standard deviation of 
the best-fit threshold fraction at a given threshold Av by substituting the relations F = bF b (Av\ At, B, a 2e — 0) and 
F = BF b into eq. with the result 

F(Av\At,a 2e =0) 

a F f it = — ; (4a ) 

y / NF(Av ta ii\At,a 2e = 0) 

This equation differs from eq.|H]in that it uses the best- fit binary fraction to infer the threshold fraction at thresholds 
At; > Avtau ■ Since the scatter in F is smaller at Av ta u , this leads to a better constraint than if the threshold fraction is 
measured directly. Comparing eqs. 1431 and IHI and using the fact that the threshold fraction F decreases monotonically 
in Av, we see that the error (Jfju < <*F as expected. Eq . 1431 approximates the error in the threshold fraction estimated 
by the Bayesian approach outlined in section [4] 
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3. APPENDIX: ANALYTIC FORMULAS FOR BINARY VELOCITY DISPERSION 

For the purpose of making analytic calculations or doing a likelihood analysis, it is useful to have a fitting function 
for the binary velocity distribution; this function was defined in section [7] (eqs. [241125]) . If the period distribution 
parameters /ii og p and <7i og p and are varied from their fiducial values, the fitting parameters remain unchanged to first 
order (although the normalizing factor changes according to eqs. |2"T1- 121)1) . Thus we have the additional benefit that 
the binary dispersion can be calculated and its dependence on the period distribution parameters is seen explicitly. In 
terms of the fitting parameters given in eqs. I24A25I and tabulated in table [1] we find the following analytic formula for 
the binary dispersion: 



Mk 2 



{/(s)erfc[-x(w ; s)\ + f(s')erk[x(v ; s')]} 



where k — 4.741 km/s is the conversion factor from AU/year, and 



/(*) 



1 



■ exp 



ClogP 



In 10 



3,s 



exp 



-2(^ log p + 3a) In 10 



x(v; s) 



3s log |v| + MiogP + 3a 



3.s 



— a logP lnlO 



(44) 



(45) 



(46) 



Although the normalization Af also varies with the period distribution parameters /ii og P and <7i og p according to 
eqs. I27ll29l this variation is weak compared to the exponential dependence in f(s) (eq. [45]). The binary dispersion en, 
is tabulated in table [1] assuming the fiducial period distribution parameters /ui og p = 2.23, <7i og p = 2.3 (P in years). 

One can also find the dispersion with a velocity cutoff v c , according to which stars with v > v c are discarded from 
the sample. We denote the dispersion with a cutoff as a(v c ). A little algebra shows that 



a 2 {v c ) = a%(v c ) + Bo%(v c ) 



where 



06 K) = r 



<tI{v c ) = 0-gerf 
v c -\v'\ 



I'c 



erfc 



V2 



<JqV2 
- erfc 



C7 V C 



V2 



vff b (log\v'\)dlog\v'\ 



(47) 
(48) 
(49) 



If a cutoff is chosen such that v c > 2.7ao, then the first term dominates in eq. |3H1 and the first term in the integrand 
of eq. |49] dominates over the second term. The complementary error function can be approximated as a step function, 
in which case the integral can be done analytically with the result 



Afn 2 



{f(s) [1 + eri\x(v ; s)\] + f(s') [etf\x(v c ; s')\ - erf\x(v ; s')\}} 



(50) 



The extra factor in eq. [50] introduced by the cutoff v c mitigates the exponential dependence of f(s) with the period 
distribution parameters (Ti og p and /iwp; indeed it nearly linearizes them over their viable range of values (roughly 
from 0-4 with P in years; see section [5f. Thus if one chooses a cutoff v c = ca where a is the measured dispersion (e.g. 
c = 2.7), the dispersion is approximately linear in the period parameters, consistent with the expected near-degeneracy 
with binary fraction. This reinforces the very useful observation that the effect of binaries on the measured dispersion 
of a galaxy can be discerned from multi-epoch data in a model-independent way. 



